PROGRAM MAPS
USE POLYMORPHIC_COMPLEXTAYLOR
TYPE(DAMAP) Id, MAP, RESIDUAL
REAL(DP) L1,KICK1,L2,KICK2,LD 
REAL(DP) :: PREC = 1.D-10
type(vecfield) F
type(ONELIEEXPONENT) SINGLE
type(DRAGTFINN) DF
type(REVERSEDRAGTFINN) RDF
INTEGER I

integer mf,me
mf=20
open(unit=mf,file='results.txt')
CALL INIT(NO1=4,ND1=2,NP1=0,NDPT1 =0)     !   <------------------ init for maps in ND1 degrees of freedom


call alloc(Id, MAP, RESIDUAL)
call alloc(F);  CALL ALLOC(SINGLE);CALL ALLOC(DF);CALL ALLOC(RDF);

KICK2=-0.9d0; L2=0.1D0


call VECFIELD_QUAD(L2,KICK2,F)

F%V(1)=0.01D0+F%V(1)
F%V(2)=0.02D0+F%V(2)

ID=1
MAP=TEXP(F,ID)

CALL PRINT(MAP,MF)

SINGLE=MAP

DO I=1,C_%ND2
 RESIDUAL%V(I)=SINGLE%VECTOR%V(I)-F%V(I)
ENDDO

WRITE(MF,*) "   "
WRITE(MF,*) " ORIGINAL VECTOR FIELD - COMPUTED ONE LIE EXPONENT "
CALL PRINT(RESIDUAL,MF,PREC)

DF=MAP

WRITE(MF,*) "   "
WRITE(MF,*) " DRAGT-FINN REPRESENTATION "
WRITE(MF,*) "   "
CALL PRINT(DF,MF,PREC)

RESIDUAL=DF
RESIDUAL=MAP-RESIDUAL

WRITE(MF,*) "   "
WRITE(MF,*) " ORIGINAL MAP - COMPUTED DRAGT-FINN  "
WRITE(MF,*) "   "
CALL PRINT(RESIDUAL,MF,PREC)


RDF=MAP

WRITE(MF,*) "   "
WRITE(MF,*) " REVERSE DRAGT-FINN REPRESENTATION "
WRITE(MF,*) "   "
CALL PRINT(RDF,MF,PREC)

RESIDUAL=RDF
RESIDUAL=MAP-RESIDUAL

WRITE(MF,*) "   "
WRITE(MF,*) " ORIGINAL MAP - COMPUTED REVERSE DRAGT-FINN  "
WRITE(MF,*) "   "
CALL PRINT(RESIDUAL,MF,PREC)


call kill(Id, MAP, RESIDUAL)
call kill(F);  CALL kill(SINGLE);CALL kill(DF);CALL kill(RDF);
close(mf);
END PROGRAM MAPS

subroutine VECFIELD_QUAD(L,K,F)
USE POLYMORPHIC_COMPLEXTAYLOR
implicit none
TYPE(VECFIELD) F
REAL(DP) L,K

! DRIFT PART
F%V(1)=L*(1.D0.MONO.'01')  /SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))
F%V(3)=L*(1.D0.MONO.'0001')/SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))

! QUADRUPOLE PART 

F%V(2)=-L*K*(1.D0.MONO.'1')
F%V(4)= L*K*(1.D0.MONO.'001')


end subroutine VECFIELD_QUAD

subroutine PBFIELD_QUAD(L,K,H)
USE POLYMORPHIC_COMPLEXTAYLOR
implicit none
TYPE(PBFIELD) H
REAL(DP) L,K

! DRIFT PART
H=L*SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))

! QUADRUPOLE PART 

H%H=H%H-(L*K/2.D0)*((1.D0.MONO.'2')-(1.D0.MONO.'002'))

end subroutine PBFIELD_QUAD

